% Figure 5. Plot the conditions for Lemma 9
% This file also computes various other comparative statics of pi
% based on US calibration
% using the model in the main text; a_L=a_H=1
% October 2022
% Makoto Nirei

clearvars;
setparams;

findq2 =@(x,rho,mu,delta,eta,pi) ...
    ((rho+mu)*delta/(eta-1))*vp(x,1-eta+mu/pi)/vp(x,mu/pi)/vp(x,1-eta)+1 ...
    -vp(x,1-eta+(rho+mu)/pi)/vp(x,-eta+(rho+mu)/pi) *vp(x,-eta)/vp(x,1-eta);    % Eq(64)

pilength=100;
piv=logspace(log10(0.002),log10(0.4),pilength);
qv=zeros(1,pilength); pbv=qv; wv=qv; thetav=qv; pstar=qv; repricefreq=qv; theta0=qv;
VL=qv; EL=qv; Tstar=qv; BB=qv; calvovar=qv; vb=qv; % preallocation
for i=1:pilength
    pi=piv(i);
    q = fzero(@(x) findq2(x,rho,mu,delta,eta,pi), [1.0000001 100]);
    pb = (vp(q,1-eta+mu/pi)/vp(q,mu/pi))^(1/(eta-1));
    w = pb*(eta-1)/eta *vp(q,eta-1-(rho+mu)/pi)/vp(q,eta-(rho+mu)/pi) * q;
    qv(i)=q;
    pbv(i)=pb;
    wv(i)=w;
    theta=vp(q,1-eta)/vp(q,1-eta+mu/pi);
    thetav(i)=theta;
    pstar(i)=q*pb;
    repricefreq(i)=mu/(1-theta) *vp(q*pb,1-eta)/vp(q,mu/pi);
    theta0_t=vp(q*pb,1-eta)/vp(q,mu/pi);
    theta0(i)=theta0_t;
    VL(i)=(theta0_t+theta0_t^2)*theta/(1-theta)^3 +theta0_t/(1-theta)^2;
    EL(i)=theta0_t/(1-theta);
    Tstar(i) = log(q)/pi;
    BB(i) = vp(q,1-eta+mu/pi)/vp(q,1-eta)/vp(q,mu/pi);
    
    % contribution of calvo on V(dlog P) (multiplied by n)
    calvovar(i) = (pstar(i)^(2*(1-eta)) - 2*pstar(i)^(1-eta) + ...
        (pstar(i)^(2*(1-eta)+mu/pi) - pb^(2*(1-eta)+mu/pi))/ ...
        ((2*(1-eta)+mu/pi) *pb^(mu/pi) * vp(q,mu/pi))) * mu/(1-eta)^2;
    
    % lower bound of v (Dec 2019)
    vb(i)=(1/rho)*(pstar(i)-wv(i))*pstar(i)^(-eta)-delta;

    % E[p^(-eta)]=N/Y
    NpYv(i)=pb*vp(q,mu/pi-eta)/vp(q,mu/pi);
    lamv(i)=mu/(q^(mu/pi)-1);
    CpYv(i)=1-delta*lamv(i);
end

% Plot Figure 5. Conditions for the upper bound of m

cap_theta=(qv.^(mu./piv)-1).*thetav;
mm=thetav.*exp(cap_theta-1)./cap_theta;

figure(7),
plot(piv,cap_theta,'-',piv,mm,'--',piv,ones(pilength,1),':','LineWidth',2)
legend('\Theta','\theta F(r*)','FontSize',16,'Location','northwest')
xlabel('\pi','FontSize',16)
ax=gca;
ax.FontSize=16;

%%------------- Dec 2021

% figure(15),
% plot(piv,wv,'b-','LineWidth',2);
% ax=gca;
% ax.FontSize=20;
% xlabel('\pi');
% ylabel('Steady-state wage rate $\bar{w}$','interpreter','latex');
% grid;
% 
% figure(16),
% plot(piv,pstar,'b-',piv,pstar./wv,'k--','LineWidth',2);
% ax=gca;
% ax.FontSize=20;
% xlabel('\pi');
% ylabel('$p^*$ and $p^*/\bar{w}$','interpreter','latex');
% legend('$p^*$', '$p^*/\bar{w}$','interpreter','latex');
% grid;
